clc;
clear;

%% main主脚本
tic;
startTime = cputime;
pde = process_pde_elliptic(1);
vector_h = [1/2 1/4 1/8 1/16];
basis_type = "P1";

u_Linf = zeros(size(vector_h,2),1);
u_L2 = zeros(size(vector_h,2),1);
u_H1 = zeros(size(vector_h,2),1);

for i = 1:size(vector_h,2)
    h = vector_h(i);
    mesh = process_mesh(pde, h);
    fs = process_fespace(mesh, basis_type);

    % process_showmesh(mesh, fs, [1 1 0]);

    A = pde.nu*term_grad_U_grad_V(mesh, fs);
    b = term_F_V(mesh, fs, pde.f);
    [A, b] = process_boundary_elliptic(A, b, mesh, fs, pde);
    U = A\b;

    process_plot_scalar(mesh, fs, U);

    u_Linf(i) = process_error_elliptic(mesh, fs, pde, U, "Linf");
    u_L2(i) = process_error_elliptic(mesh, fs, pde, U, "L2");
    u_H1(i) = process_error_elliptic(mesh, fs, pde, U, "H1");

    fprintf("已完成\t%d/%d\n", i, size(vector_h,2));
end

process_rate(u_Linf, vector_h);
process_rate(u_L2, vector_h);
process_rate(u_H1, vector_h);

TOCtime = toc;
CPUtime = cputime - startTime;
disp(["TOC time(s)", num2str(TOCtime)]);
disp(["CPU time(s)", num2str(CPUtime)]);

%% process_pde_elliptic
function pde = process_pde_elliptic(index)
switch index
    case 1
        pde.index = index;
        pde.nu = 1;
        pde.u = @(x,y) x.*y.*(1-0.5*x).*(1-y).*exp(x+y);
        pde.udx = @(x,y) 0.5*(y.*exp(y+x).*(x.^2-2.0).*(y-1.0));
        pde.udy = @(x,y) 0.5*(x.*exp(y+x).*(x-2.0).*(y+y.^2-1.0));
        pde.f = @(x,y) -y.*(1-y).*(1-x-0.5*(x.^2)).*exp(x+y)-x.*(1-0.5*x).*(-3*y-y.^2).*exp(x+y);
        pde.domain = [-1, 1, -1, 1];
        pde.bdary = [0, 1, 0, 0];
    case 2
        pde.index = index;
        pde.nu = 1;
        pde.u = @(x,y) exp(x+y);
        pde.udx = @(x,y) exp(x+y);
        pde.udy = @(x,y) exp(x+y);
        pde.f = @(x,y) -2*exp(x+y);
        pde.domain = [-1, 1, -1, 1];
        pde.bdary = [0, 0, 0, 0];
    case 3
        pde.index = index;
        pde.nu = 1;
        pde.u = @(x,y) cos(pi*x).*cos(pi*y);
        pde.udx = @(x,y) -pi*sin(pi*x).*cos(pi*y);
        pde.udy = @(x,y) -pi*cos(pi*x).*sin(pi*y);
        pde.f = @(x,y) 2*pi^2*cos(pi*x).*cos(pi*y);
        pde.domain = [0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    otherwise
        error("Invalid pde index.");
end
fprintf("pde index: %d\n", pde.index);
fprintf("pde domain: xmin,xmax,ymin,ymax\n");
disp(pde.domain);
fprintf("pde boundary: bottom,right,top,left (0=Dirichlet,1=Neumann)\n");
disp(pde.bdary);
if isfield(pde, "bdary_p")
    fprintf("pde boundary of p: (0=Dirichlet,1=Neumann)\n");
    disp(pde.bdary_p);
end
end

%% term_grad_U_grad_V
function matrix_A = term_grad_U_grad_V(mesh, fespace)
Nb = fespace.Nb;
Nlb = fespace.Nlb;
rows = zeros(mesh.Ne*Nlb*Nlb,1);
cols = rows;
val = rows;
flag = 1;
for e = 1:mesh.Ne
    V = mesh.node(mesh.elem(e,:),:);
    [int_points_loc, int_weights_loc] = integral_loc(V, "2d");
    phix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace.basis_type, "dx");
    phiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace.basis_type, "dy");
    for i = 1:Nlb
        for j = 1:Nlb
            f = phix(:,j).*phix(:,i) + phiy(:,j).*phiy(:,i);
            r = int_weights_loc * f;
            rows(flag) = fespace.elem(e,i);
            cols(flag) = fespace.elem(e,j);
            val(flag) = r;
            flag = flag + 1;
        end
    end
end
matrix_A = sparse(rows, cols, val, Nb, Nb);
end

%% term_F_V
function vector_b = term_F_V(mesh, fespace, func_f)
Nb = fespace.Nb;
Nlb = fespace.Nlb;
rows = zeros(mesh.Ne*Nlb,1);
val = rows;
flag = 1;
for e = 1:mesh.Ne
    V = mesh.node(mesh.elem(e,:),:);
    [int_points_loc, int_weights_loc] = integral_loc(V, "2d");
    phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace.basis_type, "x");
    ff = func_f(int_points_loc(:,1), int_points_loc(:,2));
    for i = 1:Nlb
        f = ff .* phi(:,i);
        r = int_weights_loc * f;
        rows(flag) = fespace.elem(e,i);
        val(flag) = r;
        flag = flag + 1;
    end
end
vector_b = sparse(rows, 1, val, Nb, 1);
end

%% process_boundary_elliptic
function [A, b] = process_boundary_elliptic(A, b, mesh, fespace, pde)
Nb = fespace.Nb;
Nlb = fespace.Nlb;
% Neumann
for k = 1:size(pde.bdary,2)
    if 1 == pde.bdary(1,k)
        bde = bdary_edge(mesh, k);
        rows = zeros(size(bde,1)*Nlb,1);
        val = rows;
        flag = 1;
        for n = 1:size(bde,1)
            points = mesh.node(bde(n,2:3),:);
            AB = (points(2,:)-points(1,:))';
            vector_n = [0, 1; -1, 0] * AB;
            vector_n = vector_n./norm(vector_n);
            fb = @(x,y) vector_n(1)*pde.udx(x,y) + vector_n(2)*pde.udy(x,y);
            e = bde(n,1);
            V = mesh.node(mesh.elem(e,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(points, "1d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace.basis_type, "x");
            ffb = fb(int_points_loc(:,1), int_points_loc(:,2));
            for i = 1:Nlb
                f = pde.nu * ffb .* phi(:,i);
                r = int_weights_loc * f;
                rows(flag) = fespace.elem(e,i);
                val(flag) = r;
                flag = flag + 1;
            end
        end
        vector_b = sparse(rows, 1, val, Nb, 1);
        b = b + vector_b;
    end
end
% Dirichlet
for k = 1:size(pde.bdary,2)
    if 0 == pde.bdary(1,k)
        bdn = bdary_node(mesh, fespace, k);
        A(bdn,:) = 0;
        A(sub2ind(size(A), bdn, bdn)) = 1;
        b(bdn) = pde.u(fespace.node(bdn,1),fespace.node(bdn,2));
    end
end
end

%% process_error_elliptic
function error_u = process_error_elliptic(mesh, fespace, pde, U, error_type)
switch error_type
    case "Linf"
        err_u = 0;
        for e = 1:mesh.Ne
            V = mesh.node(mesh.elem(e,:),:);
            [int_points_loc, ~] = integral_loc(V, "2d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace.basis_type, "x");
            uh = phi * U(fespace.elem(e,:));
            f = max(abs(pde.u(int_points_loc(:,1),int_points_loc(:,2)) - uh));
            if f > err_u
                err_u = f;
            end
        end
        error_u = err_u;
    case "L2"
        err_u = 0;
        for e = 1:mesh.Ne
            V = mesh.node(mesh.elem(e,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(V, "2d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace.basis_type, "x");
            uh = phi * U(fespace.elem(e,:));
            f = (pde.u(int_points_loc(:,1),int_points_loc(:,2)) - uh).^2;
            r = int_weights_loc * f;
            err_u = err_u + r;
        end
        error_u = sqrt(err_u);
    case "H1"
        err_ux = 0;
        err_uy = 0;
        for e = 1:mesh.Ne
            V = mesh.node(mesh.elem(e,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(V, "2d");
            phix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace.basis_type, "dx");
            phiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace.basis_type, "dy");
            uhx = phix * U(fespace.elem(e,:));
            uhy = phiy * U(fespace.elem(e,:));
            fx = (pde.udx(int_points_loc(:,1),int_points_loc(:,2)) - uhx).^2;
            fy = (pde.udy(int_points_loc(:,1),int_points_loc(:,2)) - uhy).^2;
            rx = int_weights_loc * fx;
            ry = int_weights_loc * fy;
            err_ux = err_ux + rx;
            err_uy = err_uy + ry;
        end
        error_u = sqrt(err_ux + err_uy);
    otherwise
        error("Invalid error type.");
end
end